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Abstract: Comparison of numerical performances of two methods for coefficient inverse 
problems is described. The first one is the classical Gel'fand-Levitan-Krein equation method, 
and the second one is the recently developed approximately globally convergent numerical 
method. This comparison is performed for both computationally simulated and experimental 
L_i data. 
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^ 1 Introduction 

. ^ The goal of this publication is to compare performances of two numerical methods listed in 

^ the title on both computationally simulated and experimental data. The reason why we do 

d not compare other methods is a practical one: such a study would require much more time 

and effort than the authors can afford to spend. We work with the time resolved data of 
scattering of the electrical wave field and these data are the same as ones in [21]. The data 
were collected in the field in the cluttered environment by the Forward Looking Radar of 
US Army Research Laboratory (ARL) [3D]. The goal of this radar is to detect and possibly 
identify shallow explosives. 



1.1 Previous publications 

In publication [21] the data of the current paper were successfully treated in the most chal- 
lenging case: the case of blind data. To do this, the newly developed the so-called approx- 
imately globally convergent numerical method (AGCM) was applied. Since these experi- 
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mental data are the one-dimensional ones, then it is possible to use them for a comparison 
of numerical performances of the classical Gelfand-Levitan-Krein integral equation method 
(GLK) and AGCM. "Blind" means that first the data were delivered to the mathematical 
group by engineers of ARL, Sullivan and Nguyen, who coauthor both this paper and [21]. 
No information about targets was given to mathematicians, except their placements either 
above or below the ground. Next, the mathematical group made computations by AGCM 
and delivered results to engineers. Finally, engineers compared computational results with 
the truth and informed mathematicians about their findings, see Table 1 and Conclusion in 
subsection 3.5. The most important observation was that in all available sets of experimental 
data the values of computed dielectric constants of targets were well within limits tabulated 
in [32l [33]. Since dielectric constants were not measured in experiments, then this answer 
was a quite satisfactory one. 

GLK was derived in 1950-ies by Gelfand, Levitan and Krein (THl 122]- Originally it 
was derived for the 1-d inverse spectral problem. Later, it was shown in [12] that an inverse 
problem for the 1-d hyperbolic equation can be reduced to the GLK. In addition to [12], 
we also use here the material of §4 of Chapter 1 of the book [31], where the derivation of 
[12] is reproduced. As to the numerical results for GLK, we refer to interesting publications 

AGCM was developed for Coefficient Inverse Problems (CIPs) for the hyperbolic equation 
c{x)utt = Am,x G W^,n = 1,2,3 with single measurement data, see, e.g. papers [3, IH1 El 
|20| [2H [251 [23] about this method; a summary of results is given in the book [6] . "Single 
measurement" means that the data are generated either by a single point source or by a 
single direction of the incident plane wave. Thus, single measurement means the minimal 
amount of available information. In particular, the performance of AGCM was verified on 3-d 
experimental data in [Sl[Bl[2n]- AH known numerical methods for n— d, n > 1 CIPs with single 
measurement data require a priori knowledge of a point in a small neighborhood of the exact 
coefficient. Unlike this, the key advantage of AGCM is that it delivers good approximations 
for exact solutions of CIPs without any advanced knowledge about small neighborhoods of 
those solutions. It is well known that this goal is an enormously challenging one. 

Hence, a rigorous definition of the "approximate global convergence" property was intro- 
duced in [6l [HI [211 EH]. In simple terms, it means that first an approximate mathematical 
model is introduced. The approximation of AGCM is a mild one: it amounts to the trun- 
cation of an asymptotic series with respect to l/s,s — > oo, where s > is the parameter 
of the Laplace transform with respect to t of the above hyperbolic PDE. This truncation is 
done only on the first iteration. Next, a numerical method is developed within the frame- 
work of the resulting mathematical model. Next, a theorem is proved, which claims that 
this method delivers a good approximation for the exact solution of that CIP without any 
a priori knowledge of a small neighborhood of that solution. The common perception of 
global convergence is that one should obtain correct solution if iterations would start from 
almost any point. That theorem claims, however, that a small neighborhood of the exact 
solution is achieved if iterations would start not from any point but rather from a function 
which can be obtained without any knowledge of that small neighborhood. That theorem was 
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confirmed in numerical studies of botli computationally simulated and experimental data. 
3-d experimental data were treated in [HI El [20]. In particular, the most difficult blind data 
case was considered in [HI 1201 121] • We refer to the paper [23] and references cited there for 
another non-local method for a coefficient inverse problem. 

1.2 A huge misfit between real and computationally simulated 
data 

It was pointed out in all publications about AGCM for experimental data [H [6l 1201 121] 
that there is a huge misfit between measured data and computationally simulated data, in 
the case of waves propagation in a non-attenuating medium. This misfit is evident from a 
simple visual comparison of simulated and experimental curves, see, e.g. Figures [3] and |4j 
This discrepancy is the main difficulty for the numerical treatment of experimental data by 
any computational method. In other words, the experimental data are lying far away from 
the range of the operator, which needs to be inverted to solve the CIP. It is well known that 
this is a compact operator. Hence, its range is very narrow and the inversion problem is 
unstable. Therefore, regularization is necessary. Still, the regularization theory gives recipes 
only in the case when the right hand side of the operator equation is not far from the range 
of that operator. However, in the case of those experimental data, the right hand side is 
actually far away from that range. 

Therefore, the crucial step of all above cited publications about experimental data, was 
the data pre-processing procedure. This procedure has extracted such a piece of the data 
from the whole data, which looked somewhat similar with the computationally simulated 
data. Still, it was impossible to estimate the distance between the extracted data and the 
range of that compact operator. Furthermore, those data pre-processing procedures cannot 
be rigorously justified neither from the Physics nor from the Mathematics standpoint. They 
were based on the intuition only. The single criterion of their success was the accuracy 
of resulting solutions of corresponding CIPs. Since accurate results were obtained for the 
bhnd real data in [6l [201 121] , then those data pre-processing procedures were unbiased. We 
describe below the data pre-processing procedure. It consists of two steps: extracting a piece 
of data from the whole data and multiplying this piece by a calibration factor. This factor 
is unknown and needs to be determined numerically. 

In section 2 we describe the experimental data we work with as well as the data pre- 
processing procedure. In section 3 we briefly describe the 1-d version of AGCM. GLK is 
described in section 4. In section 5 we compare performances of these two methods on our 
experimental data. We discuss results in section 6. 
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2 Experimental Data 



2.1 Data collection 

The time resolved nanosecond electric pulses are emitted by two sources installed on the radar 
(see Fig. [T|. The schematic diagram of the data collection by the Forward Looking Radar is 
depicted on Fig. [2] Only one component of the electric field is both originated and measured 
in the backscattering regime. The data are collected by sixteen (16) detectors with the time 
step of 0.133 nanosecond. Only shallow targets placed either above or a few centimeters 
below the ground can be detected by this radar. This is sufficient for antipersonnel plastic 
land mines as well as for home made explosives. The depth of the upper part of the surface of 
a shallow underground target is a few centimeters. The Ground Positioning System (GPS) 
provides the distance between the radar and a point on the ground located above that target. 
The error in the latter is a few centimeters. Time resolved voltages of backreflected signals 
are integrated over radar/target distances between 20 meters and 8 meters, and they are 
also averaged with respect to both source positions as well as with respect to readings of 
sixteen detectors. 

Since the radar/target distance is provided by GPS with a good accuracy, geometrical 
parameters of targets, including their depths, are not of an interest here. The main goal of 
the publication [24J was to calculate ratios R of dielectric constants 

^ ^ erjtarget) 
Eribckgr) ' 

where er{target) is the dielectric constant of the background medium. If er{bckgr) is known, 
then ([T]) enables one to calculate €r{target). If a target is located above the ground, then 
Eribckgr) = er{air) = 1. Unfortunately, dielectric constants were not measured during the 
data collection process. Therefore, the authors of |25 had no choice but to rely on tables 
of dielectric constants [321 [33] when evaluating the accuracy of results. In our mathematical 
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Figure 2: The schematic diagram of the data collection by the Forward Looking Radar 

model (|4])-(j8]) of section 3 i? is the function of the spatial variable x, i.e. R = R(x) : = 
G (0, 1) . Hence, after computing R (x) we took the value R to compute the dielectric 
constant of the target er{target), where 

— _ r max[o,i] if R{x) > 1, Vx G (0,1), 

~ \ min[o,i] if R{x) < 1, Wx e (0,1). ^ ^ 



Hence, in ([Ij) we actually have R := R, and by ([2j) 

Eritarget) = Reripckgr). (3) 
The importance of estimating dielectric constants of potential explosives is described in 



the following citation from the paper [21] "The recovered dielectric constant by itself is not 
a sufficient information to distinguish one target from another. The purpose of estimating 
the dielectric constant is to provide one extra piece of information about the target. Indeed, 
up to this point, most of the radar community relies solely on the intensity of the radar 
image for doing detection and discrimination of targets. It is hoped therefore that when 
the intensity information is coupled with the new dielectric information, algorithms could 
then be designed that will ultimately provide better performance in terms of probability of 
detection and false alarm rate" . 



2.2 Severely under-determined data 

The experimental data of both |24j and this paper are severely under-determined. We now 
describe main factors of the under-determination. 

Since targets can be mixtures of constituent materials and since we had worked with 1-d 
data, whereas targets are certainly 3-d, then the calculated eritarget) was a certain weighted 
average of the spatially distributed dielectric constant of a target. Due to some technical 
limitations, we possess only five sets of experimental data for five different targets. These 
targets are schematically depicted on the left panel of Fig. |3| For any target of interest, 
only a single time dependent curve can be extracted from the vast amount of data. The 
right panel of Fig. [3] depicts corresponding experimental curves, where the horizontal axis 
is time in nanoseconds. An important additional factor of the under-determination is the 
clutter surrounding targets. Other factors contributing to the under-determination were: the 
integration of all data as above, the unclear source position in ([T]), and it was also unclear 
where the time moment t = is on the data. 
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2.3 Data pre-processing 

We now describe the data pre-processing procedure. It consists of two stages. 
2.3.1 First stage: selection of a single peak 

Consider the mathematical model of the above process. Since only one component of the 
electric field was both sent into the medium and measured, and since we have received 
only one experimentally measured curve per target, we had no choice but to ignore the 
complete Maxwell's system and model the process by only one 1-d wave-like PDE. At the 
same time, recent computational testing of the full time dependent Maxwell's system in 
[TT] has demonstrated that the component of the electric field E = [Ei, E2, E^) , which was 
originally sent into the medium, substantially dominates two other components. In addition, 
we refer to the accurate performance of AGCM for experimental data in [HI El EO], where 
only a single PDE, which was a 3-d analog of (|6]), was used for the mathematical model. 

Let u{x,t), x G M, t > be the component of the electric field which is incident on the 
medium and which is also measured in the backscattering regime. Let er{x) be the spatially 
distributed dielectric constant of the medium. We assume that 

Srix) E [do,di), erEC\R), (4) 

Srix) = 1, X^(0,1), (5) 

where the numbers d^, di > 0. Thus, the interval (0, 1) is our domain of interest in the GIF. 
We consider this interval as the one whose length is one meter. We model the process of 
waves propagation via the following hyperbolic Cauchy problem 

er{x)utt = u^x, xeR, tG(0, 00), (6) 
u{x,0) = 0, Ut{x,0) = 6{x — xq), (7) 
xo = -I. (8) 

Since the source position xq is unclear in our experiment (subsection 2.2), we choose xq = — 1 
in ^ for the sake of definiteness only. 

Coefficient Inverse Problem 1 (CIPl). Determine the coefficient er{x), assuming 
that the following function g(t) is known 

uiO,t)=g{t), tG(0,oo). (9) 

The function g{t) models the backscattering data measured by the Forward Looking 
Radar. We now solve the forward problem ([gJ-Q for the case when the graph of the function 
Er (x) is the one depicted on Fig. |4]-a). The resulting function gi{t) := g(t) — H (t — \xq\) /2 
is depicted on Fig. |4]-b). Here H (t) is the Heaviside function and H {t — \x — xo\) /2 is the 
solution of the problem (6)-([8]) for the case (^) ^ 1. 

Comparison of Fig. 4 -b) with the data depicted on the right panel of Fig. |3] confirms 
the above statement about a huge discrepancy between experimental and computationally 



6 



Metal CMinder 



Taigel 1* flush Luuied metal cylnicliical cojitauiei 



Wood Stake 



Racial liiie ot siglil 




Target is wooden stake iii the gioujid 

Wood stake dielectric cojistaiit is aiomid 2 oi" 3 



Flush Buried Metal Box 

Radai' Ihie of siglil 



Taiget is fliisJi tiniietl metal b 



Bush (Clullei) 



Radar lijie of sidit 




Target is desert bush seen dicing the road 



Plastic C\ linder 



Radai hne of siglit 



Target is flush biuied empty (ah-filled) plastic 
cyhndiical contamer 

Taget-to-back-giound dielectic contiast is aiomid 1 '3=0.333 






2 6e+04 




























-2.6e+04-- 


A/ 










A// 


















ns 

























1 "tH-[)5- 
5 3e+04 




























1 








ns 


-fi.3e+()4- 





















































1 2e+04-- 
















-T.3e+03 1 












1 
































r 


s 















Figure 3: Targets and collected data before pre-processing 
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Figure 4: a) The graph of the function er{x) for which the data are computationally simulated 
via solution of the forward problem b) The computed function gi{t) := g{t) — H{t — 

\x — xq\)/2 at the edge {x = 0}, where the backscattering data are collected. The function 
g{t) is defined in ^ 



simulated data. Therefore, we must change somehow the experimental data, i.e. pre-process 
it. So that the pre-processes data would look somewhat similar with the curve of Fig. |4]-b). 
Observe that this curve has basically only one downward looking dent. If the step function 
of Fig. |4]-a) would be downward looking, then we would have an upward looking step-like 
function on the corresponding analog of Fig. |4]-b). Now, if a target is standing in the air, then 
the background medium is air with Sr^air) = 1 and er{target) > 1. Hence, it was decided 
in [21] to select only the first peak with the largest amplitude on each experimental curve. 
The rest of the curve was set to zero. Still, in order to take into account the information 
whether the target was in the air or buried, a little bit different procedure of selecting peaks 
was implemented. More precisely, the selected peak must be the earliest one with the largest 
amplitude 

all peaks for a buried target, 
all downward looking peaks for a target in the air. 

Also, since the time moment {t = 0} was unknown to us, we set t := to be such a point 
on the time axis, which is one (1) nanosecond to the left from the beginning of the selected 
peak. Resulting superimposed pre-processed curves are shown on Fig. |5} It is clear now 
that the condition t E (0, oo) in ^ is not a restrictive one, because each curve of Fig. |5] 
represents a function with a compact support for t G (0, oo) . It is also clear from comparison 
of Figures |4]-b) and [5] that the pre-processed data are far from the range of the operator 
A {sr) '■= u (0, t) , which should be inverted to solve CIPl. 
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Figure 5: Pre-processed experimental data from the right panel of Fig. 3. 
2.3.2 Second stage: the choice of the calibration factor 

The modulus of the curve of Fig. |4]-b) does not exceed 0.2. On the other hand, amplitudes 
of peaks of experimentally measured curves of Fig. [3]are of the order of 10^. This means that 
the data correspond to A5 {x — xq) in ([T]), where the number A = const. > is unknown. In 
our experience a similar situation takes place for all experimental data, not necessarily for 
the current one. It is clear therefore that in order to make the data suitable for the model 
(|6|-([9]), they must be multiplied by a certain number, which we call "calibration factor" and 
denote CF. A proper choice of CF is a non-trivial task. To have unbiased studies, it is 
important that CF must be the same for all targets. 

In [24j the same calibration number CF was chosen for all five targets. It is important 
that because of the blind case, that choice was unbiased. If we would have measured the 
dielectric constant of at least one target in advance, then we would choose CF in such a way 
that the number i? in ([T]), which would be obtained via solving the above CIP, would be close 
to the measured one. However, since dielectric constants were unknown in advance, then 
another method of choice of CF was used in [23]. Namely, consider the Laplace transform 
of the function u {x, t) , 

oo 

=£(„,(x,«)):=/„,(..<)e-'*. s = const. >0. (10) 



We call the parameter s pseudo frequency. Denote 

9 is) = Cg, (11) 
f \ ^ riu(+ I exp(-s|xo|) 

9{s) = g{s)-go{s). 
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Let 'gsim (s) be the Laplace transform of the function depicted on Fig. |4]-b). Let 'gbush (s) 
be the Laplace transform of the function, which corresponds to bush on Fig. [5j Then we 
multiply all experimental data by such a calibration number CF that values of functions 
^sim (s) and CF ■ 'gbush (s) would be not far from each other for s G [1,5]. It was shown 
in [2lj that CF = 10~^ was good for this goal. Furthermore, we have discovered that 
the same value of CF = 10~^ provides similar behavior for Laplace transforms of four 
other curves of Fig. [s] In [23] we have also varied the calibration factor CF by 20% as 
CFi = 0.8 ■ 10-^ CF2 = 1.2 ■ 10"^ In all three cases of CF, CFi and CF2 the resuhing 
values of dielectric constants were within tabulated limits [23] ■ Hence, we have chosen in 

m 

CF = 10-\ (12) 
The choice of CF for the GLK is a more delicate issue, and we discuss it in section 5.1. 

3 Approximately Globally Convergent Method 

In this section we briefly describe the 1-d version of AGCM for CIPl. Details can be found 
in [23] as well as in section 6.9 of the book [6j. It was proved in [2^ that the function w{x, s) 



in (10) is the solution of the following problem 

Wxx — s'^Er (x) w = — 5 (x — Xo) , x G M, Vs > s, (13) 
lim w{x, s) = 0. (14) 

|a;|— >oo 

Denote 

w{0,s) ■.= g{s) , s>s, (15) 



where the function g (s) is defined in (11). It was shown in [23] that 



Wx{0, s) := p{s) = sg (s) - exp(sxo). (16) 



3.1 Integral differential equation 



Consider the fundamental solution wo{x,s) of the problem (13), (14) for er{x) = 1. Then 



2s 

Since w{x, s) > [23], we can consider In [w{x, s)] . Denote 



Wo{x,s) = — exp(— s|x — Xol). (17) 



r(x,.) = ^"[^^/"^;^^^'^^] . (1^ 



Then (13)-(17) imply that 

Txx + s^rl - 2srx = e^(x) - 1, x > 0, (19) 
r(0,s) = if.is), r,(0,s) = ¥Pi(s), (20) 



10 



f \ ^ n -t \ wo M , ^0 f \ 2 exp(sxo) 
fois) = In^(s) - ln(2s) + — , <^^{s) = - 



r,(l,s) = 0. 



Differentiate equation (19) with respect to s. Denote 

q{x,s) = dsr{x,s), ^/'o(s) = (^^(s), tp^is) = ip[{s), 

s 

r{x,s) = — J q{x,T)dr + V{x,s), 

s 

V{x,s) = 



{In [w{x, s)] — In [wo{x, s)]} = r{x, s) 



(21) 
(22) 

(23) 
(24) 

(25) 



Here s > is a sufficiently large number which is chosen in numerical experiments. We call 
V{x, s) the tail function. Actually s is one of regularization parameters of this method. We 
have 



V{x,s) = + O ( ^ ) , s CX), 



q{x,s) 



s 

p{x 



(26) 
(27) 



Using (19)- (25), we obtain 



s s 

q^^-2s'^q^ J q^{x,T)dT + 2s jq 



-I 2 



(x, r)(ir 



s 

2sq^ + 2 J q^{x,T) dr 



s 

+ 2s^q^V^ - 4s V; j q^ix, T)dr + 2sV^ - 2V; = 0, s G [s, s], 



(28) 



q{Q,s) = iI)q{s), g^(0,s) = V'i(s), g^(l,s) = 0, se[s, s]. (29) 
The main difficulty of this numerical method is in the solution of the problem (p8|, 



(|29|). Equation (|28|) has two unknown functions q and V. To approximate both of them, 

. First, given an 



we use a scheme, which is based on both inner and outer iterations 
approximation for V , we update q via solution of a boundary value problem on the interval 
X G (0, 1). This is our inner iteration. Next, we update the unknown coefficient er{x) and 



solve the forward problem (13), (14) for s := s on the entire real line M with this updated 



coefficient er{x). Next, we update the tail function V {x, s) via (25). If both functions q and 



V are approximated well by functions q and V respectively via this numerical procedure, 
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then we use backwards calculations (24), (19) to approximate the unknown coefficient er{x), 
r{x,s) = — J q{x,T)dT + V{x,s), (30) 

s 

~ , , f l + r^x + s^7l-2s?^, a;e(0,l), 



3.2 Numerical procedure 

Consider a partition of the interval [s, s] into small subintervals with the grid step size 
/i > and assume that the function q (x, s) is piecewise constant with respect to s, 

s = Sn< Sn-1 < ■■■ < So = s, Si-i - Si = h] q {x, s) = g„ (x) , for s G (s„, . 

For each subinterval (s„, we obtain a differential equation for the function g„ (x) . We 
assign for convenience of notations go — 0. For each n we iterate to improve tails. Next, we 
go to the next n. This way we obtain functions g„,fc, Ki,fc- The equation for the pair (g„,fc, Ki,fc) 
is 



n-l 



n" 



'n-l 



j=0 
2 



n-l 



n-l 



(32) 



J2q'A +2hJ2q'^+2A,,^V;JhJ2l- 
0=0 / j=0 \ j=0 

-A2,n{V;,k)' + '^A2,nV;^k, X G (0, 1) , n G [1 , iV] , A; G [1 , m] , 
where the number m should be chosen in numerical experiments. The boundary conditions 



are generated by ( 29 ) 



qn,k (0) = ^Ao.n, q'n^k (0) = ^i,„, qik (1) = 0. 



(33) 



Boundary conditions (33) are over-determined ones. On the other hand, our attempts to use 



only two out of three boundary conditions (33) never led to acceptable results [23]. Thus, 



given the function ^ (x) , we use all three conditions (33) and find an approximate solution 



(ln,k (x) of the problem (32), (33) via the Quasi- Reversibility Method (QRM), see subsection 



3.3. We refer to the book [26j for the originating work about the QRM and to [HK 
for some follow up publications. QRM is well suited for solving over-determined boundary 
value problems. Convergence of QRM is proved via Carleman estimates |21| . 

The choice of the first tail function Vi^i (x) is described in subsection 3.3. Let n > 1. 
Suppose that for j = 0, ...tt. — 1 functions qj (x) , Vj (x) are constructed. We now construct 
functions qnk,Vnk for k = l,...,m. We set Vni{x) := Vn-i{x). If the function Vnk{x) is 
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given, then we find an approximate solution qn,k (x) of tlie problem (32), (33). Next, we find 
the approximation e^"''^^ for the unknown coefficient Er (x) via discrete analogs of (30), (31), 



n-l 



rn,k (x) = -hqn,k -h'^qj + Vn,k, x E (0, 1) , 



j=0 



1 + <,k {^) + 4 Kk i^)] ' - 2s„r;^, {x),xe (0, 1) 
l,a;^ (0,1). 



Next, we solve the forward problem (13), (14) with Er (x) := Er 



Sn,k) 



X] 



s and find 



the function Wn,k+i {x,s) this way. After this we update the tail via formula (25) in which 
w{x,s) := Wn,k+i {x,s) . This way we obtain a new tail Vn,k+i (x) . Similarly we continue 
iterating with respect to tails m times. Next, we set 



Qn (x) 



(x) , K (x) := (x) := K,™ (x) , 4") (x) := 4"'"^) (x) 



replace n with n + 1 and repeat this process. By [23J we continue this process until 



either 



^2(0,1) 



<10-^or ||VJ„(g„,fc)||^.o^,)>10 



(34) 



where the functional Ja {qn,k) is defined in (35) of subsection 3.3. Here the norm in the space 



L2 (0, 1) is understood in the discrete sense. In the case when the second inequality (34) 
is satisfied, we stop at the previous iteration taking e^"''^ (x) as our solution. If neither 
of two conditions (30) is not reached at n := iV, then we repeat the above sweep over the 



interval [s, s] taking the tail Vn (x) as the new tail Vi^i (x) . Usually at least one of conditions 



(34) is reached at most on the third or on the fourth sweep, and the process stops then. 



3.3 Computing functions g„ fc(a:), Vi 

As it was mentioned in Subsection 3.2, we compute approximations for functions qn.k (x) via 
QRM. Let L {qn,k) (x) be the operator in the left hand side of equation (32) and P„^fc (x) its 



right hand side. Let a G (0, 1) be the regularization parameter. The QRM minimizes the 
following Tikhonov regularization functional 



\Ln,k ijln,l 



1 1 2 I II II 2 



(35) 



subject to boundary conditions (33) 



We now describe an important step of choosing the first tail function Vi^i (x) . Consider 



the asymptotic behavior (26) and (27) of functions V {x,s) and g(x,s) with respect to the 
truncation pseudo frequency s — )■ 00. We truncate in (26) and (27) terms 0(l/s^) and 
0(1/ s^) respectively. This is somewhat similar with the definition of geometrical optics 
as a high frequency approximation of the solution of the Helmholtz equation. Hence, for 
sufficiently large s 



V (x, s) 



p [Xj 



q{x,s) 



p (x) 



(36) 
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Hence, setting in equation (|28j) s 
equation for the function p (x) 



s and using (36), we obtain the following approximate 



p" (x) = 0, X e (0, 1) . 



(37) 



Boundary conditions for the function p (x) can be easily derived from (29) and (36) as 

p (0) = -iV (s) , p' (0) = -s^V-i (s) , p' (1) = 0. 



(38) 



We find an approximate solution p (x) of the problem (37), (38) via the QRM. Next, we set 
for the first tail function 

p (x) 



Via ix) 



(39) 



Theorem 1 is a simplified version of Theorem 6.1 of [23], also see Theorem 6.7 of ^ for 
the 3-d case. 

Theorem 1. Let the function e* (x) satisfying conditions Q), he the exact solution 
of our CIP for the noiseless data g* (t) in Fix the truncation pseudo frequency s > 1. 
Let the first tail function Vi^i (x) he defined via (31)-(39). Let 7 G (0, 1) he the level of the 
error in the houndary data, i. e. 

IV^o is) - i^l {s) I < 7, (s) - ^* (s) I < 7, for s G [s, s] , 



where functions ^/'q (s) , t/']^ (s) depend on the function g{t) in ^ via (11), (21), (23) and 
functions ip^ (s) , ipl (s) depend on the noiseless data g* (t) in the same way. Let ^Ja = 7 
and h = max(7, /i). Let Q he the total numher of functions e-r^'''^ computed in the ahove 



algorithm. Then there exists a constant D = 
is so small that 

h < 



D (xo, do, di, s) > 1 such that if the numher h 

(40) 



1 



then the following convergence estimate is valid 



I (n,fc) _ * I 



(41) 



where the numher u G (0,1) is independent on n,k,h,e^(^''^\e*. 

Theorem 1 guarantees that if the total number Q of computed functions er"''^^ is fixed and 
error parameters cr, h are sufficiently small, then iterative solutions e(^'^'^ (x) are sufficiently 
close to the exact solution e*, and this closeness is defined by the error parameters. Estimate 



(41) guarantees the stability of AGCM with respect to a small error in the data. Therefore 



the total number of iterations Q can be considered as one of regularization parameters of our 
process. Two other regularization parameters are the numbers s and a. The combination 



of inequalities (40) and (41) has a direct analog in the inequality of Lemma 6.2 on page 



156 of the book [15] for classical Landweber iterations, which are defined for a substantially 
different ill-posed problem. Indeed, it is stated on page 157 of the book [15j that the number 
of iterations can serve as a regularization parameter for an ill-posed problem. 
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The main advantage of Theorem 1 is that even though it does not require any advanced 
knowledge of a small neighborhood of the exact solution e*, it still guarantees the most 
important property of this algorithm: the property of delivering at least one point in this 



neighborhood. In principle, (|41|) guarantees that one can take the function resulting 



from just the first tail Vi.i (x) as an approximate solution. In other words, already on the first 
iteration we obtain a good approximation. Next, one can proceed either of two ways. First, 
one can continue iterations as long as the total number Q of iterations satisfies inequality 



(40). In doing so, one can choose an optimal number of iterations, based on numerical 
experience, see (34) for our stopping rule. Second, one can choose either of functions ei'^''^-' 



with nk < Q as an approximate solution and then refine it using one of locally convergent 
numerical methods. This is a two-stage numerical procedure, which was developed in [HI [6] 
and some other publications about AGCM. An effective method of this sort is, e.g. Adaptive 
Finite Element Method (adaptivity), which minimizes the Tikhonov functional on a sequence 
of locally refined finite element meshes, see, e.g. |5| [10] and Chapter 4 of the book fB]. 



In the numerical implementation of AGCM the functional {qn,k) in (35) was written 
via finite differences and was minimized with respect to the values of the function g„ ^ at 
grid points. Conjugate gradient method was used for the minimization. Convergence of this 
method to the unique minimizer of Jq, {qn,k) follows from the strong convexity of Jq, {qn,k) , 
which, in turn follows from Lemma 5.2 of f23]. The s— interval was s G [s,s] = [0.5, 12] . The 



length of each small subinterval (s„, s„_i) was h = 0.5. We took a = 0.04 in (35). For each 
n = 1, ...,N we have calculated functions g„,A: for k = 1, 10 := m. We refer to [231 El] for 
details of the numerical implementation of AGCM for this case. 

3.4 Numerical results for AGCM 

Test 1. We now test AGCM for computationally simulated data for the case of the function 
€r (x) displayed on Fig. |4]-a). We have introduced the multiplicative random noise in the 
function g (t) in ^ as 

g^iU)=giU)il + aQ, (42) 

where t = tj is the discrete value number i of the variable t G (0, T), which was used in data 
simulations, a G (0, 1) is the noise level and is the value number i of the random variable 
^ G (—1, 1). Hence, e.g. a = 0.05 means 5% of the noise level in the data. Because of Fig. 
|4]-b) as well as because of the rapid decay of the kernel e~'^* of the Laplace transform (10) 



we took T = 4. Images with 0% , 5% and 10% noise in (42) are displayed on Fig. [6] One can 
observe that both the location of the abnormality and the maximal value of the coefficient 
Er (x) are imaged with a good accuracy for all three cases. 

Test 2. We now show the image using experimental data for case of "bush". This case 
is the most difficult one, since bush is highly heterogeneous. We remind that because our 
data are severely under-determined ones, the calculated function R (x) in ([T]) is a sort of a 
weighted average value of the ratio ([TJ over the volume of a target. Next, we have applied 
formulas (|2]), ([3]) to find €r (target) . We do not reproduce images of other four available cases 
of real data, since they are published in [21]. Fig. M displays this image. Hence, R = 6.4 
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(a) 



0,2 0.4 




(b) 




(c) 



0,6 0,8 



Figure 6: Images from computationally simulated data for different values of the random noise 
in the data, see (42). The true target is the same as on Fig. ^a). a) 0% noise, b) 5% noise, c) 
10% noise. 



7r Er 




0.0 0.25 0.5 0.75 1.0 



Figure 7: The blindly computed function R{x) for bush, see the fourth row of Fig. Q The 
computed value maxi?(x) = R = 6.4 = ^^(bush) is within tabulated limits of [3, 20] /74f 

in this case. As one can see from Fig. [3| bush was standing in the air. Hence, in this 
case Er ihckgr) = 1. We conclude therefore that the computed value Er (bush) = 6.4. On the 
other hand it is clear from [H] that the tabulated value is Er (bush) G [3, 20] . Therefore, our 
blindly computed value of the dielectric constant of bush is within tabulated limits. 

Table 1 [23] shows computed values for all five cases of experimental data. The back- 
ground value of the dielectric constant of the dry sand, in which three out of five targets 
were buried, was between 3 and 5 [32] • It is clear from Fig. [3] which targets were buried 
and which ones were standing in the air. The dielectric constant for metals is not defined 
in Physics. Nevertheless, it was established in numerical experiments in [21] that as soon 
as the value of the dielectric constant of a target Er > 10, the backscattering electric wave 
from this target is about the same as the one from a metallic target. Therefore, the so-called 
appearing dielectric constant for metals was introduced in [23], i-e. this is such a dielectric 
constant, which is between 10 and 30. Thus we assign for this constant 

Er (metals) G [10, 30] . (43) 
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Tabulated values of dielectric constants of the wood stake and plastic cylinder in Table 1 
are taken from Tables [321 [33]. As to the bush, see [H]. The tabulated value of (metal) is 
taken from (43). Note that published values of dielectric constants are given within certain 
intervals. Therefore, by ([3]), values of computed and tabulated (target) in Table 1 are also 
given within certain intervals denoted by brackets. This is unlike computed values of R. 



Table 1: Blindly computed dielectric constants for five targets of Fig. [3] as in (12), by the 
globally convergent numerical method versus tabulated ones. "1" in the second column 
means that the target was in the air and "[3,5]" means that it was buried in the dry sand 
whose dielectric constant is between 3 and 5 1321. 



target 


R 


Srihckgr) 


computed Sr (target) 


tabulated Ej. (target) 


metal cylinder 


4.3 


[3,5] 


[12.9,21.4] 


[10,30], see ( 


43 


) 


metal box 


3.8 


[3,5] 


[11.4,19] 


[10,30], see ( 


43 


) 


wood stake 


3.8 


1 


3.8 


[2,6], see [3^ 




bush (clutter) 


6.4 


1 


6.4 


[3,20], see [H] 


plastic cylinder 


0.28 


[3,5] 


[0.84, 1.4] 


1.2, see [33j 



Conclusion. Since dielectric constants of targets were not measured in experiments, the 
only thing we can do is to compare computed values of Sr (target) with tabulated ones. The 
most important conclusion, which can be drawn from Table 1, is that in five out of five blind 
cases computed values of dielectric constants of targets were well within tabulated limits. 



We remind that the variation of the calibration factor CF =10 ^ in (12) by 20% has also 
resulted in values being well within tabulated limits. 

4 The Gelfand-Levitan-Krein Integral Equation (GLK) 

4.1 Forward and inverse problems 

The forward problem in this case is 

er{x)utt = u^x, X > 0, t > 0, (44) 
u{x,0) =ut{x,0) = 0, (45) 
u.\x=o = 5it). (46) 



Coefficient Inverse Problem 2 (CIP2). Determine the function er{x) in (44), assuming 
that the following function f{t) is known 

u{0,t) = f{t), t>0. (47) 
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Without loss of generality, we assume that 

£,(0) = 1. 

Change the spatial variable as 



X Z = r(x) = J a/ Er {T)dT. 




(4J 



(49) 



Hence, taking into account (|48|), we obtain that (|47|) holds at 2; = and (44)- (46) imply that 

(50) 
(51) 



Utt = u^z —Uz, z > 0, t> 0, 



a[z) 

u {z, 0) = Uz {z, 0) = 0, 



Uz\z=0 

a{z) 



er{r-\z)) 



(52) 
(53) 



Consider even extensions of functions a{z) and u{z^ t) in {z < 0}. To do this, we set without 
changing notations a{z) = a{—z),u{z,t) = u{—z,t) for z < 0. It follows from [12] that we 
obtain then the following Cauchy problem 



a'iz 



Utt 



(T{Z 



-Uz, z G 



t > 0, 



Hence, 
where the function uiz, t 



m(2,0) = 0, Ut{z,Q) = -25{z). 

u{z,t) = -^/a{z) ■ Hit - \z\) + u{z,t), 



(54) 



(48), we obtain 



is continuous for {z,t) G M x (0,T) ,VT > and u{z,t) = for 

cr (0). Combining this with 

(55) 



t < 1^1 . Hence, (|53|) and pAj imply that /(0+) = u(0+, 0' 

/(o+) = -i. 



Let 



fit) 



fit), t>0 
-fi-t) t<0' 

Hence, the function fit) is odd and its first derivative /'(t) is even. Choose a number T > 0. 



Taking into account (55) and using results of [T2l ITS] , we obtain that the GLK equation has 
the form 



wiz,t)-^ I f'it-r)wiz,T)dr=^, te[-z,z], WzE[0,T/2]. 



(56) 
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This is not a Volterra equation, because z G [0,T/2] is a parameter here. Hence, equation 



(56) should be solved for all values of the parameter z G [0, T/2]. The solution of this equation 
is connected with the function cr(2;) via 



lim w {z, t) :- 



ze[0,T/2]. 



Hence, by (53) and (57) 



(r-' (-)) 



16w^{z,z ) 



(57) 



(5^ 



Therefore, in order to solve CIP2, one needs first to solve equation (56) for all z G [0,T/2]. 
Next, one needs to reconstruct the function Sr (t~^ (z)) via (58). Finally, one needs to change 



variables backwards transforming z into x via inverting (49). 

Theorem 2 [121 EI]- Let the function f{t) be given for t G [0,T] for a certain value 
of T > 0. Suppose that there exists unique function er{x) of CIP2 satisfying conditions 
Q), in which M is replaced with {x > 0} . Also, assume that condition (48) holds. Let 
t{z) be the function defined in (49). Then solution of the GLK integral equation (56) exists 
for each z G [0,T/2] and is unique. And vice versa: if the function f{t) is known for all 
t G [0, T] and the integral equation has unique solution for each z G [0,T/2], then CIP2 
has the unique solution er{x) for x G [0,r^^ (^/2)] satisfying above conditions. 



4.2 Performance of GLK for computationally simulated data 



In this subsection we present results of our computations for GLK for computationally 
simulated data. For our testing we choose the coefficient Sr (x) whose graph is depicted 
on Fig. |4]-a). This is the same function as the one for AGCM. We use the FDM to solve the 
forward problem (44)- (46). The solution of this problem generates the data f{t) := u{0,t), 
t G [0, T], where T > is a certain number. To verify the stability of our computations with 



respect to the noise in the data, we introduce the random noise in the data as in (42), where 
g it) is replaced with / (t) and T is replaced with T. 



Next, we use the change of variables x -v^ 2; in (49) and define the time interval t G 
(0,T),T = T{T). Indeed, since in (56) z G [0,T/2] and t G [— 2;, z] , then the value of 
-2max = T/2 should be such that (^max) := Xmax > 1- Note that all what we need is 
a lower estimate of 2;max rather than its exact value. Hence, to get such an estimate, it is 
sufficient to know a lower estimate for the unknown coefficient Sr (x) rather than the function 
Sr (x) itself. Let 

{Zn}n=0 C [0. ^/2]> = Za<Zr<...<ZM^= T/2 

be a partition of the interval z G [0,T/2] into small subintervals with the grid step 
size hz = Zn — Zn-i- For each Zn we approximate the integral (56) via the trapezoidal rule. 
For the variable t G [—Zn, Zn] we use discrete values {tk = ±-2^}^""^ . Next, we solve the 
resulting linear algebraic system for each Zn and find the discrete function w{zn,tk),Wk G 



-n,n] ,Vn G [0,M^ 



Next, we find the discrete function (t ^ i^n)) via (58). Finally, we 
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for a certain partition {a;„}„; 



numerically invert ( 49 ) and find the discrete function e,. 
of the interval [0,Xinax ^ [0, 1] • This discrete function Sr (a;„) is an approximate solution of 
CIP2. 

Since the differentiation of noisy data is an unstable procedure and f'{t — A) is the kernel 
of the integral equation (56), we calculate the derivative f'{t) in a special way. For each 
point tk ■= Zk,k = 1 



Mz — 1 we approximate the derivative f'{tk) as 



fit 



k+l) 



2K 



(59) 



M- 



Next, we extend this approximation as the discrete even function for = ~Zk, k = 1, 
1. Even though (59) seems to be a non-regularizing procedure, it works quite well for our 



goal, since in the discrete integration in (56) the right hand side of (59) is actually multiplied 
by K/2. 



We test GLK for 0%, 5% and 10% noise in the data f(t). Let and be the number 
of grid points and the mesh size respectively for the variable x G (0, 1). We took = 250, 
which led to = 0.004. After the change of variables (49) we got z G [0, 1.168], M^, 



_ 349, 
3.34 ■ lO""*. Results of computations are displayed on Fig. |8[a)-f). Note that Figures 
c), e) depict the function 



fit) 



4 



(fit) + 1) 



(60) 



Comparison of Figures [sj-b), d), f) with Figures [6]-a), b), c) shows that GLK better 
images the shape of the target than AGCM. The inclusion/background contrast maxe.,. (x) 
is computed accurately by both GLK and AGCM. On the other hand, the contrast computed 
by AGCM does not change when the level of noise changes, while it does change for GLK. 



5 Comparison of AGCM and GLK on Real Data 

As it was pointed out in subsection 2.3.2, to work with experimental data, it is necessary 
to figure out the calibration factor. This factor was established for AGCM as CF = 10~^, 



see (12). Since dielectric constants of targets were not measured in experiments, while 
CF = 10~^ still led to acceptable values of those constants, then we believe that the best is 
to use computed values of R in Table 1 to figure out the calibration factor for GLK. Thus, 
we proceed as described in Subsection 5.1. 



5.1 The choice of the cahbration factor for GLK 

The procedure of this choice consists of the following three steps: 

Step 1. Select any target X listed in Table 1. We call X "calibration target". Let Fx (t) 
be the data for the target X, which are pre-processed as in subsection 2.3.1, see Fig. |5j 

Step 2. Let Rglk [X] and Ragcm (X) be the function R (x) in ([T]) computed by GLK and 
AGCM respectively for the target X. Let Rglk i^) and Ragcm be corresponding values 
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(I.I) 0.2 4 06 0.8 1.0 



Figure 8: Performance of GLK for the computationally simulated data for the target of Fig. ^a ). 
a), c) and e) the computed function f{t) in (60) for 0%, 5% and 10% noise respectively, b), d) 



and f) computed function er{x) for cases a), c) and e) respectively. 



of the number R in (|2j). Hence, Ragcm (X) is the number R for the target X in Table 1. 
Let CF (X) be the cahbration factor for GLK. Thus, GLK should be computed for the data 
CF (X) ■ Fx (t) . Choose such a calibration factor CF (X) that Rglk {X) = Ragcm (X) . 
Thus, Rglk {X) is computed via GLK for the data CF (X) ■ Fx (t) and Ragcm {X) is 



computed via AGCM for the data 10""^ Fx (t) , see Q. 

Step 3. Let F 7^ X be any other target listed in Table 1. When applying GLK to the 
data Fy (t), use the same calibration factor CF (X) . 
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Figure 9: a) The function Rglk{x) for the wood stake when it was the calibration target. 
Here CF(wood stake) = 10~^, Rglk = 3.8; b) the function i?GLir(plastic cylinder) for the 
case when the plastic cylinder was the calibration target. Here CF(plastic cylinder) = 1.8 ■ 10~^, 
i?GLx(plastic cylinder) = 0.28; c) the function RcLKi^^^^) for the case when bush was the 
calibration target. Here CF(bush) = 0.6 ■ 10~^ ^cLi^ (bush) = 6.4. 

5.2 Results 

When applying GLK to the pre-processed data, we have chosen T to be the moment of time 
which is 1.33 nanoseconds to the right from the end of the selected peak. 

Case 1. Wood stake is the cahbration target. Using Step 2 of Subsection 5.1, we have 
obtained CF(wood stake) = 10~^. Figure [9[a displays the function Rglk (wood stake). 



Table 2: Case 1. Computed numbers Rglk for four other targets in comparison with the 
values of R taken from Table 1. Here wood stake was the calibration target. Bold faced is 
an unrealistic value of Rglk- 



Target 


Rglk 


R from Table 1 


metal cylinder 


3.25 


4.3 


metal box 


3.97 


3.8 


plastic cylinder 


0.6 


0.28 


bush 


42.5 


6.4 



We observe from Table 2 that the value Rglk (bush) = 42.5 is significantly overestimated 
compared with the published interval R G [3, 20] in |14j . 

Case 2. Plastic cylinder is the calibration target. We have obtained CF(plastic cylinder) = 
1.8-10^^. Figure[9]-b displays the function i?Gix(plastic cylinder) with this calibration factor. 
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Table 3: Case 2. Computed numbers Rglk for four other targets in comparison with values 
of R taken from Table 1. Here plastic cylinder was the calibration target. Bold faced are 
those values of Rglk which result in unrealistic values of dielectric constants. 



Target 


Rglk 


R from Table 1 


metal cylinder 


9.8 


4.3 


metal box 


13.9 


3.8 


wood stake 


17 


3.8 


bush 


»100 


6.4 



We observe from Table 3 that the values of Rglk for the wood stake and bush are 
significantly overestimated, compared with the published intervals [2, 6] of [32] and [3, 20] of 
[H] respectively. Actually we had a blow in the image of bush via GLK in this case. 

Case 3. Bush is the calibration target. We have obtained CF(bush) = 0.6 ■ 10~^. 
Figure |9]-c) displays the function i?G'Li^(bush) with this calibration factor. 



Table 4: Case 3. Computed numbers Rglk for four other targets in comparison with values 
of R taken from Table 1. Here bush was the calibration target. Bold faced are those values 
of Rglk which lead to unrealistic values of dielectric constants. 



Target 


Rglk 


R from Table 1 


metal cylinder 


2.0 


4.3 


metal box 


1.16 


3.8 


wood stake 


2.2 


3.8 


plastic cylinder 


0.74 


0.28 



It follows from Table 4 that the value of Rglk for metal box is significantly underes- 
timated. Indeed, since in the dry sand, which was the background here, Sr G [3,5] [32] . 
then, using (|2]), we obtain (metal box) G [3.48,5.8]. On the other hand, by (43) we 
should have Er (metal box) G [10, 30]. Thus, we have the underestimation of Sr (metal box) 
by the factor of 1.72. In addition, the value of Rglk for plastic cylinder is overestimated 
in Table 4. Indeed, since (background) G [3,5], then it follows from this table that 
(plastic cylinder) G [2.22,3.7]. On the other hand. Table [33] tells us that we should 
have Er (plastic cylinder) ^ 1.2. Hence, we have overestimation of Er (plastic cylinder) by the 
factor of 1.85. 

Our calculations for the case when metal box was chosen as calibration target led to the 
same results as those in Table 2. This is because Ragcm (wood stake) = Ragcm (metal box) = 
3.8. Also, Ragcm (metal cylinder) = 4.3 ~ Ragcm (metal box) . Hence, the calculation for 
the case when the metal cylinder was chosen as calibration target again led to results similar 
with those of Table 2. 
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Figure 10: Dependencies of Ragcm 3nd Rglk from the calibration factor for bush and wood 
stake are displayed. In both a), b) "1" for AGCM corresponds to CF = 10^" , see ( [T^ . a) "1" 
corresponds CF(bush) = 0.6 ■ 10"^ (Case 3). b) "1" corresponds CF(wood stake) = 10"^. 



To see how the changes in the cahbration factor affect the reconstructed values of Ragcm 
and Rglk, we present Figures 10 The visual analysis of these curves shows that the increase 
of the calibration factor affects Ragcm linearly, and it affects Rglk exponentially. This 
explains results of above Cases 1-3. 



6 Discussion 

Theorem 2 guarantees existence and uniqueness of the solution of the GLK equation only 
when the data f{t) belongs to the range of the operator of the forward problem, i.e. for the 
case of errorless data. However, in the realistic case of an error in the data, stability with 
respect to this error is not guaranteed by Theorem 2. Hence, the question about regularizing 



properties of GLK remains open. On the other hand, estimate (41) of Theorem 1 ensures 
stability of AGCM with respect to a small error in the data, and the same is true for the 
3-d version of this method [HI El 125] . 

We have used both computationally simulated and experimental data to compare nu- 
merical performances of two methods: GLK and AGCM. In the case of computationally 
simulated data, the above difference between Theorems 1 and 2 is slightly reflected in the 
different behavior of computed inclusion/background contrasts in the presence of noise: these 
contrasts were not changing with noise in the data for AGCM and were slightly changing 
for GLK. On the other hand, in the case of the above synthetic data, GLK provides better 
images of shapes of targets, compared with AGCM. 

As it was pointed out in subsection 2.3.2, to have unbiased studies of these experimental 
data, it is important that the calibration factor should be chosen the same for all five targets 
of Table 1. However, it was shown in subsection 5.2 that it is impossible to choose such a 
calibration factor for GLK, which would provide satisfactory values of inclusion/background 
contrasts for all five targets. On the other hand, AGCM has worked successfully with the 
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uniform calibration factor CF = 10~^ for all five targets in the case of blind data. We 
point out that these conclusions are relevant only for the specific case of the above set of 
experimental data and for the above specific data pre-preprocessing procedure. On the other 
hand, we do not possess other sets of 1-d experimental data at the moment, and we are also 
unaware about other appropriate data pre-processing procedures. 
To explain Fig. 



10, we present the GLK equation (56) as 



w{z, t) 



CF (X) 



f'{t — r)w{z, T)dT 



1 

2' 



t G [-z,z], \/z G [0,r/2]. (61) 



In (61 ) f{t) is the pre-preprocessed data as in in subsection 2.3.1 and f{t) is the odd extension 



of f(t). Denote f3 := z ■ CF (X) sup|i|<j^/2 /'(^) • If £ (0, 1), then one can solve integral 



equation (61) via the resolvent series, see, e.g. the book |35]. It is clear from this series that 



the solution w{z,t) changes almost with an exponential speed when the calibration factor 
CF {X) changes. This explains, at least partially, the exponential behavior of GLK curves 



of Fig. 10 This is because the data for the inverse problem are introduced in the operator. 



On the other hand, almost linear behavior of AGCM curves on Fig. [10] can be explained 

see 



by the fact that AGCM uses the logarithm of the solution w of the problem (13), (14) 



(18). 
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